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Metastability in a nano-bridge based hysteretic DC-SQUID embedded in 
superconducting microwave resonator 

Eran Segevjj Oren Suchoi, Oleg Shtempluck, Fei Xue|l| and Eyal Buks 

Department of electrical engineering, Technion, Haifa 32000, Israel 

(Dated: July 30, 2010) 

We study the metastable response of a highly hysteretic DC-SQUID made of a Niobium loop 
interrupted by two nano-bridges. We excite the SQUID with an alternating current and with direct 
magnetic flux, and find different stability zones forming diamond-like structures in the measured 
voltage across the SQUID. When such a SQUID is embedded in a transmission line resonator similar 
diamond structures are observed in the reflection pattern of the resonator. We have calculated the 
DC-SQUID stability diagram in the plane of the exciting control parameters, both analytically 
and numerically. In addition, we have obtained numerical simulations of the SQUID equations of 
motion, taking into account temperature variations and non-sinusoidal current-phase relation of the 
nano-bridges. Good agreement is found between experimental and theoretical results. 

PACS numbers: 85.25.Dq, 74.40.Gh, 74.25.Sv 
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I. INTRODUCTION 

In the last two decades Superconducting Quantum In- 
terference Devices (SQUIDs) have regained the inter- 
est of researchers worldwide, due to their use as solid- 
state quantum bits. More recently, such SQUIDs were 
embedded in superconducting Transmission Line Res- 
onators (TLRs), in order to produce circuit cavity quan- 
tum electrodynamics in the strong [l|-|6| and the disper- 
sive [3-[2| coupling regimes [10|. Other applications, in 
which SQUIDs are mainly used as nonlinear classical el- 
ements, include Josephson bifurcation amplifiers [lll-[l^ 
and tunable resonators [15l - [l7j . Tunable resonators were 
also used to demonstrate parametric amplification and 
squeezing [lo-^l] , and might also be used to demonstrate 
the dynamical Casimir effect [22, 23]. 

The injected power into a TLR is usually limited by 
the critical current of the DC-SQUID embedded in the 
resonator. The upper bound of this current is determined 
by the sum of the two critical currents of the Josephson 
junctions (JJs) composing the DC-SQUID. While typi- 
cal critical currents of DC-SQUIDs range between few to 
tenth of microamperes, the applications that involve res- 
onance tuning and parametric amplification could benefit 
from larger critical currents. In order to have larger crit- 
ical currents one can fabricate DC-SQUIDs using nano- 
bridges instead of JJs. Nano-bridges, which are merely 
artificial weak-links having sub-micron size, were shown 
to have similar Current-Phase Relationship (CPR) as JJs 
under certain conditions [24-26]. These nano-bridge JJs 
(NBJJs) are characterized by l arge critical current, on 
the order of milliamperes [27',^. DC-SQUIDs with large 
critical currents are often characterized by hysteretic re- 
sponse and metastable dynamics |29l-l33|. These char- 
acters are naturally made extreme in NBJJ based DC- 
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SQUIDs. In addition, NBJJs have very high plasma fre- 
quency [34], on the order of one Terahertz, which enables 
operation of the microwave TLR without introducing in- 
terstate transitions in the embedded SQUID. Thus one 
could employ a SQUID as a nonlinear lumped inductor, 
operating at microwave frequencies. 

In this paper we experimentally and numerically study 
metastable response of NBJJs based DC-SQUID, sub- 
jected to an alternating biasing current. We first the- 
oretically analyze stability zones of a highly hysteretic 
DC-SQUID in the plane of the bias current and mag- 
netic flux control parameters. Then we directly measure 
the voltage across a SQUID in that plane. Comparison 
between experimental results and between analytical and 
numerical theoretical predictions yields good agreement. 
Moreover, we measure the reflection spectra from several 
devices integrating a SQUID and a TLR and find qualita- 
tive agreement between the theory and the experiments. 



II. EXPERIMENTAL SETUP 

A simplified circuit layout of our devices is illustrated 
in Fig. [IJb). We fabricate our devices on high resistiv- 
ity Silicon wafers, each covered by a thin layer of Sil- 
icon Nitride. Each device is made of Niobium, having 
layer thickness of less than 100 nm, and is composed of 
a stripline resonator having a DC-SQUID embedded in 
its structure. The resonator is designed to operate in the 
gigahertz range, having a length of about / = 19 mm, 
which sets its first resonance mode around 2.5 GHz. The 
DC-SQUID has two NBJJs, one NBJJ in each of its 
two arms (see Fig. [21(b) and inset of Fig. [3|). The 
NBJJs are fabricated using FEI Strata 400 focused ion 
beam system [27, 35] at an accelerating voltage of 30 kV 
and Gallium ions current of 1.5pA. The outer dimen- 
sions of the bridges range from 100 x 100 nm^ for large 
junctions to 60 x 80 nm^ for relatively small junctions. 
The actual dimensions of the weak-links are smaller be- 
cause the bombarding Gallium ions penetrate into the 
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FIG. 1: (Color online) (a) Measurement setup, (b) Schematic 
(unsealed) layout of our device. Red and blue background col- 
ors indicate experimental setup having temperatures of 300 K 
and 4.2 K, respectively. 



Parameter 


E19 


E38 


E42 


SQUID Type 


RF 


DC 


DC 


SQUID Area [/um^] 


1936 


870 


1057 


Nb Thickness [nm] 


50 


100 


60 


Self-Inductance [pH] 


127 


112 


141 


/c[mARMs] 


— 


2.29 


1.74 


/3l Calc 


— 


1089 


106 


/3l Fit 


— 


722 


83 


^LFit 


1.5 


— 


45,35 


a Fit 


— 


0.026 


0.032 



TABLE I: SQUIDs parameters. The self-inductance was nu- 
merically calculated using FastHenry computer program 39| . 
The parameter /3l Calc was evaluated analytically using the 
measured critical current. The parameters /3l Fit, 13 l Fit, 
and a Fit were evaluated according to fittings of stability di- 
agrams to measured data. 



Niobium layer, and consequently suppress super conduc- 
tivity over a depth estimated between 30 nm and 50 nm 
123, S, [33- Despite the small dimensions of our NBJJs, 
most of our SQUIDs have critical currents on the order of 
milliamperes (see table H]). A feed-line, weakly coupled to 
the resonator, is employed for delivering input and out- 
put signals. An on-chip transmission line passes near the 
DC-SQUID and is used to apply magnetic flux through 
the DC-SQUID at frequencies ranging from DC to the 
gigahertz. An on-chip filtered DC bias line is connected 
directly to the DC-SQUID and is used for direct mea- 
surements of the SQUID. The Low-Pass Filters (LPFs) 
are designed to minimize the degrading effect of these 
connections on the quality factor of the resonator. Some 
measurements are carried out while the device is fully 
immersed in liquid Helium, while others are carried out 
in a dilution refrigerator where the device is in vacuum. 
Further design considerations and fabrication details can 
be found elsewhere [34, 38]. 

The experimental results in this paper are obtained 
from three devices whose parameters are summarized in 
table m The experiments are carried out using the setup 
depicted in Fig. [TJa). We report on two types of experi- 
ments. In the first one we obtain low frequency current- 
voltage measurements of the SQUID using the DC bias 
line, while the resonator does not play any role. We 
use a lock-in amplifier, which applies alternating current 
through the SQUID, having excitation frequencies on the 
order of kilohertz. We measure the voltage across the 
SQUID using the lock-in amplifier and, in addition, we 
record the spectral density of the voltage using a spec- 
trum analyzer and its time domain dynamics using an 
oscilloscope. In the second type of experiments we inves- 



tigate the response of an integrated SQUID-TLR device 
to a monochromatic incident probe tone that drives one 
of the resonance modes of the TLR. The reflected power 
spectrum is recorded by a spectrum analyzer. In such ex- 
periments the DC bias line is left floating, and does not 
play any role in the measurement. In both types of exper- 
iments we apply DC magnetic flux through the SQUID, 
and in experiments with TLRs we also add modulated 
magnetic flux at gigahertz frequencies. 



Numerical method 

Simulations of DC-SQUID circuit model (see Fig. 
[21(a)) are done by numerically integrating its Equations 
Of Motion (EOMs); Eqs. © and (|8|). We introduce 
a sinusoidal excitation current to the EOMs and calcu- 
late the phases of the two NBJJs composing the DC- 
SQUID, the DC-SQUID voltage versus time, and the 
Fourier transform of this voltage at the frequency of ex- 
citation. The experimental measurement frequency was 
usually in the order of 1 kHz, which is about nine orders 
of magnitude smaller than the SQUID plasma frequency, 
and about six orders of magnitude smaller than the time- 
scale of thermal processes in the SQUID [40]. Therefore, 
in order to make the simulations feasible in terms of com- 
putation time, we have made two ease assumptions. The 
first assumption is that the excitation frequency used in 
simulation can be increased as long as the DC-SQUID 
follows this excitation adiabatically. Adiabatic approx- 
imation of NBJJs based DC-SQUID is thoroughly ana- 
lyzed in Ref. [34], where it is shown that the plasma 
frequency of a NBJJs based DC-SQUID is expected to 
be in the order of 1 THz. Thus in practice, the excita- 



(a) 



i__flfYYn 

, — nnnm 



80nm 



" LI 2 



R,%Ic,%^C, 



^j> hi 



^ w 






o. 



Q 



u 1/2 

L- nnrnn- 




FIG. 2: (a) Circuit model of a DC-SQUID, (b) Scanning elec- 
tron microscope (SEM) image of a nano-bridge, whose dimen- 
sions are 80 x 60 x 50 nm^ . 



I2 are the currents flowing in the upper and lower arms 
respectively (Fig. Efa)). The second is external mag- 
netic flux <l>x applied through the DC-SQUID. The total 
magnetic flux threading the DC-SQUID loop is given by 
<l> = $x + !//_, where /_ = (/i — I2) /2 is circulating cur- 
rent in the loop. Assuming sinusoidal CPR, the Joseph- 
son current /j/e in each junction {k = 1,2) is related to 
the critical current I^k and to the Josephson phase 7^ 
by the Josephson equation /j^ = I^k^^'^^k' By employ- 
ing the coordinate transformation 7+ = (71 + 72) /2 and 
7_ = (71 — 72) /2, and the notation I^ = {Id + /c2) and 
/c- = (^ci — ^02), one finds that the potenti al g overning 
the dynamics of the DC-SQUID is given by |31|] 



■ cos 7+ cos 7- 



a sm 7+ sm 7_ 



tion frequency used in simulation is set between 100 MHz 
and 1 GHz. The second assumption is that thermal pro- 
cesses have a negligible influence on the experimentally 
measured dynamics of the SQUID, as long as the excita- 
tion frequency is kept low. Therefore, although we use 
high excitation frequency in the simulations according to 
the first assumption, the temperature of the SQUID is 
held at base temperature and does not evolve with the 
SQUID dynamics. On the other hand, in simulations 
related to measurements done with an integrated TLR- 
SQUID device, is which the measurement frequency is 
high, thermal effects are taken into account by including 
thermal EOMs; Eq. 1^. 



III. THEORY OF HYSTERETIC DC-SQUID 

In this section, we develop a theory of hysteretic and 
metastable DC-SQUID, taking into account the self- 
inductance and asymmetry of the DC-SQUID. Similar 
models were also developed by others, but us ually lack- 
ing comparison with experimental data [2_9|, l32|, l4lM45| or 
the emphasis put on other aspects of SQUID dynamics 
[3l|, |46|, |47|. We begin with a lossless model, in order to 
extract the SQUID stability diagram, and add the dis- 
sipation and fluctuation terms to the EOMs on a later 
stage. 

The circuit model is shown in Fig. [2fa) (a SEM image 
can be seen in the inset of Fig. [3]). It contains a DC- 
SQUID having two NBJJs, one in each of its arms, with 
critical currents Id and /c2, which in general may differ 
one from another. Both NBJJs are assumed to have the 
same shunt resistance Rj and capacitance Cj (which is 
considered extremely small for NBJJs [48]). The self- 
inductance of DC-SQUID, L, is assumed to be equally 
divided between its two arms. Typical inductance values 
of our SQUIDs are listed in table HI 

The DC-SQUID is controlled by two external parame- 
ters. The first is bias current /x = ^1 + ^2, where /i and 



+ 7- + 



^0 



/Pl 



-7+ 



(1) 



where /3l = ttLIc/^o is a dimensionless parameter char- 
acterizing the DC-SQUID hysteresis, a = Ic-/Ic char- 
acterizes the DC-SQUID asymmetry, £^0 = (^o^c) /27r is 
the Josephson energy, and <l>o is flux quantum. 



A. Stability Zones 

The extrema points of the DC-SQUID potential are 
found by solving 



du 



57+ 
du 

97- 



/x 

sm 7+ cos 7_ + Of cos 7+ sm 7_ — — = 


: , (2) 


27_ 
cos 7+ sm 7_ + a sm 7-1- cos 7_ + —-— 

PL 


27r$, 



0. 



(3) 



In general, Eqs. (J2j) and (J3j) have periodic solutions, 
where the solutions differ one from another by 27r7Ti+ in 
7+, 27rm_ in 7_, and by 2$o^- in ^x, where m+ and 
rri- are integers. 

The Jacobian of the potential u is given by 



J = 






c^7_ c^7-)- c^7^ 



(4) 



For local minima points of the potential, both eigenvalues 
of J are positive. Thus, we find boundaries of stability 
regions of these minima points in the plane of the JJ 
phases, 7+ and 7_, by demanding that 



det J = 0, tr J > 0. 



(5) 



Furthermore, the matrix J is independent of both con- 
trol parameters I^ and <l>x, thus also these boundaries are 
independent of I^ and $x- Finding the stability thresh- 




FIG. 3: (Color online) Potential diagram of a DC SQUID, 
i^(7+,7_) , drawn using /x = 0.1/c, and $x = 0. Local min- 
ima points are labeled by black dots and saddle points are 
labeled by plus signs. The inset shows a SEM image of a 
DC-SQUID. 



olds in the plane of I^ and ^x is done by substituting the 
solutions of Eq. ^ in Eqs. (J2]) and (J3]). 

Figure [3] plots the potential diagram of a DC-SQUID, 
calculated for E42 using /3l = 83 and a = 0.032. The 
solution of Eq. (|5]) produces the black closed curves that 
enclose the DC-SQUID local minima points in the plane 
of 7+ and 7_. The local minima points are labeled by 
black dots and saddle points are labeled by plus signs. A 
local minimum point losses its stability when increasing 
either |/x| or \^^\ to a point where it merges with one of 
the saddle points close to it. 



The limits of small and large /3-l 

In general, solving Eqs. (J2]) and (J3]) for a given set 
of 7+ and 7_ can only be done numerically. Analytical 
solution can be derived only for the extreme limits of 
/3l ^ 1 and /3l ^ 1. In the former limit, which is not 
the focus of this paper, the derivation leads to the well 
known formula for the SQUID critical current [49| 



1 - (1 - a2)sin^ 



^0 



Ic 



(6) 



In the opposite limit of /3l ^ 1, the condition of det J = 
implies that cos 71 cos 72 = 0, and the condition of 
tr J > implies that (1 + a) cos 71 + (1 — a) cos 72 > 0. 
Consider first the solution for the minimum point near 
(7+,7_) = (0,0). Other solutions can be obtained from 
the periodic properties of Eqs. (J2j) and (J3j). For this so- 
lution, Eq. [5] is satisfied along a square that is formed 
by the lines connecting the four vertexes 71 = ±7r/2 and 



72 = ±7r/2. Substituting these vertexes into Eqs. (J2j) and 
([3]) yields a bounding contour which has a rectangle shape 
with vertexes at (/x/^, 27r<l>x/^o/^L) = (1,<^), (o^,!), 
(—1, —a) and (—a, —1) in the plane of the control param- 
eters /x and <l>x (see Fig. lU^a)). This rectangle crosses 
the axis /x = at the points <l>x = ±<I>o/3l (1 — o^) /27r 
and the axis $x = at the points /x = ±/c (1 — o^). 



B. Equations of motion 

Applying Kirchhoff's laws on the DC-SQUID cir- 
cuit model, substituting Josephson's current-phase and 
voltage-phase equations, and taking into account the fluc- 
tuation dissipation theory, yields the following EOMs for 
the SQUID phases 71 and 72 [49| 

71 + /^d7i + (1 + «o) y (6i) sin 71 

+ ^^ (71 - 72 + 27rc&x/^o) = /x//cO + ^nl, (7) 
PLO 



72 + /^d72 + (1 - Q^o) y (©2) sin 72 

-^^ (71 - 72 + 27r^x/<^0) = /x//cO + ^n2, (8) 

PLO 

where the overdot denotes a derivative with respect to 
a normalized time parameter r = cjpit, where uo^\ is the 
SQUID plasma frequency, and /3d = 1/ (i^jCjcjpi) is the 
damping coefficient. In general, the NBJJ critical cur- 
rent, /c/c (^ = 1,2), and thus I^ and /c-, are tempera- 
ture dependant, as we discuss later. Thus we employ the 
notation /co/c, ^co, and /co- for the corresponding critical 
currents at a base temperature Tq, which is the tempera- 
ture of the coolant. In addition, we employ the notation 
/3lo = ttL/co/^o and a^ = (/coi + ^02) / (^oi - ^02). 
The term y{Qk)^ where O^ = T^/Tc is the normal- 
ized temperature of the k^^ NBJJ, expresses the depen- 
dence of the NBJJ critical current on its temperature, 
and is equal to unity as long as the temperature of the 
NBJJs are held at base temperature. The factor g^k is 
a noise term, whose spectral density for the case where 
hv/k^T < 1 is given by Si^ {v) = 4.kBT/Rj, with A:b be- 
ing the Boltzmann constant. In what follows we neglect 
this noise term in the numerical simulations. 

To evaluate the voltage across the DC-SQUID, denoted 
as l^sQD, we assume that the loop inductance is equally 
divided between its two arms, and get 



VsQD = 



1 


"$o 


/d7i d72^ 
{dt dt ^ 


\ Ldl^ 
1 2 dt 


2 


Un 



(9) 



It is known that NBJJs may have complex CPR 
[25|, l50|, Isjj, which deviates from the normal sinusoidal 
CPR of a regular J J. According to a theory, brought 
in the appendix for completeness, such deviation would 
modify the sin 7^ term in the SQUID EOMs, Eqs. ([7]) 
and ([8|), which become Eqs. ()A5p and ()A6p . We have 
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FIG. 4: (Color online) (a) Stability diagram in the plane of 
the control parameters /x and $x, drawn for /Sl = 20 and 
a = 0.1. The curves are drawn using solid line for /x > and 
dashed line for /x < 0. The number close to each curve la- 
bels the number of flux quanta trapped in the DC-SQUID in 
the corresponding LSZ. Corresponding numbers are also indi- 
cated in panel (c). The bold black curve marks the threshold 
between static and oscillatory zones. The double-headed blue 
arrow is drawn between the points /x = =b0.9/c, $x = 0.52$o, 
and marks the control parameter range used in the numerical 
simulation shown in panels (b) and (c). These panels show 
SQUID voltage in the time domain (b) and the NBJJs phases 
in the phase space (c). The time scale is normalized by the 
period of excitation Tx = cJx/Stt. 



made some simulations in which moderate changes in the 
CPR of our SQUIDs is assumed, and found no significant 
difference between these results and results that neglect 
this deviation. In addition, the physical dimensions of 
our NBJJs are relatively small, and the measured values 
of /3l of our SQUIDs are relatively high, thus following 
the explanations detailed in appendix Al of Ref. [34], the 
effect of non-sinusoidal CPR in our NBJJs is expected to 
be negligibly small. Therefore, our conclusion is that our 
SQUIDs can be modeled by normal sinusoidal CPR. 



C. Stability Diagram 

The stability diagrams in the plane of the control pa- 
rameters /x and <l>x, and in the plane of the DC-SQUID 
phases 7+ and 7_, are plotted in Fig. HJ^a) andJU^c), re- 
spectively. The curves were first calculated in the plane 
of 7+ and 7_ by numerically solving Eq. (J5j) , with the pa- 
rameters /3l = 20 and a = 0.1. The above solutions were 
then substituted into Eqs. (J2j) and ([3]) to produce the 
closed contours of panel (a). Each closed contour bounds 
a Local Stability Zone (LSZ) corresponding to a different 
integer number of fiux quanta trapped in the DC-SQUID. 
Corresponding LSZs in Fig. |4] (a) and HJ^c) are labeled 
by the same number, indicating the number of trapped 



fiux quanta. 

The stability diagram can be separated into two global 
zones. The first is called static zone [46|, |47|, where the 
system has one or more LSZs depending of the value of 
the hysteresis parameter /3l. The static zone is bounded 
in the horizontal axis by a threshold current called the 
oscillatory threshold, given by /th (*^x) = h — I- (^x)- 
This threshold is periodic on the external fiux <l>x, hav- 
ing a maximum value equal to the critical current of the 
DC-SQUID. When the DC-SQUID is biased to the static 
zone, it is always found in a LSZ, though transitions be- 
tween LSZs may be forced by the control parameters. 
The second zone, called oscillatory zone [46|, [43 (also 
known as dissipative or free-running zone), spreads over 
two unbounded areas for which the excitation current is 
larger (in absolute value) than the oscillatory threshold. 
In this zone the DC-SQUID has no stable state; it oscil- 
lates at very high frequencies and dissipates energy. In 
what follows we focus our study on the static zone. Fur- 
ther study of the dynamics in the region of spontaneous 
oscillations can be found in Ref. 1471. 



D. SQUID Dynamics 

The EOMs, Eqs. ([7|) and (|8|), were numerically in- 
tegrated using the parameters which used to draw the 
stability diagram of Fig. IH and using the control param- 
eters <l>x = 0.52^0 and I^/Ic = 0.9. The range of the 
excitation is marked by a double-headed, blue arrow in 
the stability diagram of Fig. IH^a). Note that the left 
arrow head crosses the threshold separating LSZ_0 and 
LSZ_1 (LSZs corresponding to or 1 trapped fiux quanta, 
respectively) ; whereas the right headed arrow crosses the 
threshold separating LSZ_1 back into LSZ_0. The results 
of a simulation in which a DC-SQUID is periodically ex- 
cited along this path are shown in Fig. HI Panel (b) 
shows the DC-SQUID voltage in the time domain. Each 
excitation cycle contains two spikes, which occur close 
to extrema points of the excitation amplitude. Panel 
(c) shows the simulation in the phase plane of 7+ and 
7_. The simulation shows that the system periodically 
switches between LSZ_0 and LSZ_1, and that the dynam- 
ics does not involve any additional LSZs. Thus a positive 
spike in the time-domain response corresponds to a tran- 
sition from LSZ_0 to LSZ_1, and a negative spike corre- 
sponds to an opposite transition. When the excitation 
is monotonically increased the system mostly lingers in 
LSZ_0 and when it is decreased it mostly lingers in LSZ_1. 
Note that a serial resistor R = 0.5 Q was post-simulation 
added to the voltage results in order to emphasize the ex- 
citation cycle. Such a serial resistance also unavoidably 
exists in the wiring of our experimental setup. Note also 
that the duration of each spike, which is related to the 
relaxation time of the SQUID, is negligible compared to 
the period of excitation, Tx = uj^f^TT. Thus, the assump- 
tion that the response of the DC-SQUID to the excitation 
is adiabatic is reasonable. 
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FIG. 5: (Color online) Experimental (a) and numerical cal- 
culations (b) of direct voltage measurements and calculations 
of E38, in the parameter plane of /x and $x- Black contours 
mark the corresponding stability diagram, drawn using solid 
line for /x > and dashed line for /x < 0. The diagram is 
applied only on half of the colormap in order to leave some 
data uncovered. The white bold contour marks the PDSZ 
threshold. Black bold contour marks the oscillatory thresh- 
old. Each plus sign marks a set of control parameters, and 
the corresponding number counts the LSZs boundaries that 
the SQUID periodically crosses if excited by that set of pa- 
rameters. 



E. Periodic dissipative static zone 

Figure [5] shows a quantitative comparison between ex- 
perimental results measured with E38 using lock-in am- 
plifier (panel (a)) and simulation results calculated using 
the corresponding parameters (panel (b)). Both panels 
show colormaps of the voltage across the DC-SQUID as 
a function of control parameters, /x, and <l>x. The black 
contours in both graphs mark the corresponding stabil- 
ity diagram. The negative part of the stability contours 
is folded onto the positive part and is drawn by dashed 
lines. Thus, solid lines represent stability thresholds for 
positive excitation currents, and dashed lines represent 
stability thresholds for negative excitation currents. The 
bold black contour is the threshold between the static 
zone and the oscillatory one. 

The folding of the stability contours divides the static 
zone of the stability diagram into regions having diamond 
shapes in the plane of I^ and <l>x. Each diamond re- 
gion bounds a range of parameters for which the system 
crosses the same number of stability thresholds during 
an excitation cycle. For example, when excited with pa- 
rameters bounded by the diamond marked by 2, the DC- 
SQUID would cross a single threshold during the positive 
duration of the excitation cycle and another one during 
the negative duration, thus returning to the original LSZ, 
similar to the case discussed in Fig. HJ^c). Each crossing 



FIG. 6: (Color online) Simulation results for E38. Panels (ai), 
where i = 8, 14 corresponding to the marked points in Fig. 
m plot the DC-SQUID voltage in the time domain, panels 
(hi) magnify the corresponding marked squares in panels (ai). 
Panels (ci) show the simulation results and the corresponding 
stability zones in the phase space of 7+ and j- . 



of a threshold line, either by positive or negative currents, 
triggers a spike, which in turn contributes additively to 
the measured voltage across the DC-SQUID. Thus, each 
diamond bounds a range of parameters for which a sim- 
ilar voltage is measured. The diamonds follow the peri- 
odicity properties dictated for the control parameters by 
Eqs. (0) and daj). 

We define an additional threshold, marked by bold 
white contour, which separates the static zone into two 
sections. In the first, which applies for periodic excita- 
tion currents smaller (in absolute value) than this thresh- 
old, the DC-SQUID is captured in a single LSZ after fi- 
nite number of excitation cycles, and no spikes in the 
voltage appear afterwards. We call this section Periodic 
Non-Dissipative Static Zone (PNDSZ). In the second sec- 
tion, which spreads for excitation currents larger than the 
white threshold but smaller than the oscillatory thresh- 
old, the SQUID periodically jumps between LSZs and 
dissipates energy. Thus, we call this section Periodic Dis- 
sipative Static Zone (PDSZ), and call the white threshold 
itself the PDSZ threshold. 

Figure [6] shows results of two simulations, calculated 
using the control parameters marked by points 8 and 14 
in Fig. [5l Panels (az) where i = 8, 14 show the DC- 
SQUID voltage as a function of time, calculated during 
two excitation cycles. Each cycle has a bunch of spikes 
during its positive duration and another bunch during 
the negative one. Panels {hi) magnify the corresponding 
bunch of spikes marked by dashed squares in panels (az). 
In panel (b8) one counts four spikes corresponding to 
four crossings of stability thresholds. Panel (c8) shows 
the phase space dynamics in the plane of 7_ and 7+. 
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FIG. 7: (Color online) Experimental measurements obtained from E42. Panels (a) and (d) draw the first and second SQUID 
voltage harmonics respectively. Panel (b) draws the voltage noise level at frequency of IkHz. Panel (c) draws a statistical 
analysis of time-domain voltage traces, which calculates the difference between probabilities of measuring positive and negative 
voltage spikes respectively. 



The system periodically cycles between five LSZs, where 
most of the time it lingers either in the upper left or in 
the lower right LSZs. The transition between these LSZs 
is forced by the driving sinusoidal bias current, which 
induces four jumps, corresponding to the four spikes in 
panel (b8). Panels (bl4) and (cl4) show a richer dy- 
namics, which emerges from the fact that point +14 is 
located inside the oscillatory zone. The dynamics in- 
cludes both forced transitions between LSZs and sponta- 
neous oscillations related to the oscillatory zone. These 
two kinds of transitions are clearly seen in Fig. [5l panel 
(cl4), which plots the dynamics of 7+ and 7_ in the 
phase space. The first kind, in which 7_ and 7+ change 
monotonically, corresponds to forced transitions between 
LSZs. The second kind, in which the system cranks be- 
tween LSZs, corresponds to the spontaneous oscillatory 
dynamics. The spontaneous oscillations last as long as 
the temporal driving bias current is greater than the os- 
cillatory threshold. The distinction between the forced 
transitions and the spontaneous oscillations can also be 
noticed in the time domain, shown in Fig. [5l panel (bl4), 
where the first six spikes are similar and distinct, whereas 
the rest of the spikes emerge in pairs, in which one spike 
is slightly stronger than the other. 



F. Hybrid Oscillatory Zones 

Figure [71 shows low-frequency experimental measure- 
ments of E42 DC-SQUID, which has a moderated self- 
inductance parameter /3l = 83 (compared to the one of 
E38). Similar to the experiment with E38, we excite 



the SQUID by a sinusoidal current having frequency of 
cjx/27r = 2.5 kHz, and measure the voltage across the 
DC-SQUID as a function of excitation current I^ and 
external magnetic flux <l>x. In addition, we measure the 
spectral density of the voltage using a spectrum analyzer 
and its response in the time domain using an oscilloscope. 
Figure [71(a) shows this voltage and has the correspond- 
ing folded stability diagram drawn on it. The colormap, 
which focuses on the oscillatory threshold, reveals an in- 
teresting phenomenon. The threshold related to positive 
excitation currents (solid bold line) is shifted compared 
to the one related to negative excitation currents (dashed 
bold line). This mismatch between the thresholds cre- 
ates hybrid oscillatory zones, in which the DC-SQUID is 
driven to the oscillatory zone either for positive currents 
or negative currents, but not for both. 

Figure [TJb) plots a colormap of the voltage noise level, 
measured using the spectrum analyzer at a frequency of 
1 kHz (any frequency which is not an integer harmonics 
of cjx gives similar results). Plotting the noise level is a 
good way to detect thresholds since noise-rise is gener- 
ally expected near any bifurcation threshold [52]. Two 
distinct patterns of noise-rise are clearly seen in the col- 
ormap. One is related to the solid hue positive threshold 
and the other to the dashed line negative threshold. Note 
that this noise-rise is almost undetectable in the lock-in 
measurements. 

The existence of hybrid zones is easily observed in time 
domain measurements. Figure [HI shows two pairs of time 
domain traces, each pair has an experimental trace (pan- 
els (le), (2e)) and a simulated trace (panels (Is), (2s)), 
measured (calculated) for the parameters labeled by plus 




FIG. 8: Time domain experimental measurements (left col- 
umn) and numerical simulation (right column) for E42. The 
two rows correspond to the two sets of control parameters 
marked in Fig. [71 



the measured voltage is asymmetric in time, i.e. nei- 
ther symmetric nor anti- symmetric time response. The 
hybrid zones are characterized by such a time response, 
and indeed, the plotted colormap shows strong response 
of the second harmonics in those zones. 



IV. PARAMETRIC EXCITATION 

In recent years several demonstrations of using 
SQUIDs to manipulate the resonance fre quen cies of a su- 
perconducting resonator were reported p^5l-[l7l. l3^ . A 
SQUID in such applications is usually considered as a 
nonlinear variable inductor embedded in the resonator 
in a way that couples the resonance frequencies of the 
resonator to the SQUID impedance. The variation of 
the SQUID inductance is usually done by changing the 
magnetic flux through the SQUID, whereas the current 
through the SQUID is defined by the state of the coupled 
system. 



A. Stability Zones 



marks and the corresponding number in Fig. [Tl^a) and 
[71(b). The first (second) pair is related to the hybrid zone 
where only positive (negative) currents drive the system 
to the oscillatory zone. Accordingly, the first (second) 
pair experiences only positive (negative) spikes. The dif- 
ference between the measured and simulated line shapes 
of the spikes might be due to finite (about MHz) band- 
width of our measurement setup which is far too low 
to resolve high frequency spikes. Therefore, measured 
spikes in the output signal merge into one continuous and 
slowly decaying pulse. Local heating effects, which are 
neglected in this simulation and will be discussed later, 
may also degrade DC-SQUID performance and suppress 
the spikes. 

A statistical analysis of the time-domain behavior is 
seen in Fig. Efc), which plots the probability of counting 
a positive spike PJ^^®, minus the probability of counting 
a negative one Pg"®^, during a single lock-in excitation 
cycle. The probabilities were measured for the param- 
eter range spanned by I^ and ^x, and were calculated 
on 2 sec long traces with cJx/Stt = 2.5 kHz. This col- 
ormap clearly reveals the existence of the hybrid oscilla- 
tory zones, where the differential probability of counting 
spike having either positive or negative polarities is al- 
most one. The direction of the spikes agrees with the pre- 
diction of the stability diagram. Outside the hybrid zones 
the colormap shows near zero differential probability. In 
these areas there are no spikes at all, or the counting 
of positive and negative spikes is similar. This behavior 
could be used for creating bidirectional DC-SQUID sen- 
sors in which the polarity of measured voltage indicates 
the polarity of the detected flux change. Panel (d) plots 
the second harmonics of the measured SQUID voltage (at 
5 kHz). This harmonics is expected to be amplified when 



In the following, we analyze the stability of a DC- 
SQUID, excited by a magnetic flux having constant and 
alternating parts. Consider the case where ^x is given 

by 



$, = $o(*r + *rcosK,i)) 



(10) 



where $^^/$o and <l>^^/<l>o are arbitrary amplitudes. Re- 
call that for the case of /3l ^ 1, and for the minimum 
point near (7^,7-) = (0,0), the bounding rectangle 
crosses the axis /x = at the points ^^ = ±<l>o/3L/27r, 
where /3l = (1 — a) /3l. The range of stability for the 
minima points near (7+, 7_) = (nvr, nvr), where n is inte- 
ger, is given by 



27r $0 27r 



(11) 



Furthermore, the stability condition is achieved when 
the largest value of ^x, i-e. ^^^ + ^J^, coincides with the 
largest value of the stability range, i.e. j3i,/2t:^ or when 
the smallest value of <l>x, i.e. ^"^ — $5^, coincides with 
the smallest value of the stability range, i.e. — /3L/27r. 
Thus the boundary contours in the plan of ^^^ and ^^ 
are given by the two equations 



^dc ^ ^a 



ZTT 



(12) 



In experiments with resonators we employ the para- 
metric amplification method of operation |18l-[21|. We 
inject a relatively weak probing signal into the resonator, 
having frequency equals to one of the resonance frequen- 



cies of the resonator, cjx, and measure the reflected power 
off the resonator using a spectrum analyzer. Note that 
the DC connections do not play any role in this measure- 
ment. A weak signal is a one for which the current gener- 
ated through the SQUID is much smaller than its critical 
current. In addition, we applied constant and variable 
magnetic flux through the SQUID, given by Eq. (p!Q|) 
with cjpx = (2cJx + ^^)^ where Auj is taken to be much 
smaller than the resonance bandwidth of the resonator. 
The measured reflected power spectrum includes a tone 
at cjx and several sidebands spaced by ±772 • Auj from cjx, 
where m is an integer. These sidebands are the products 
of the nonlinear frequency mixing between uj^ and ujp^. 
Although this mixing process is more complex than our 
direct SQUID measurements, it should essentially follow 
the same SQUID dynamics, provided that adiabatic ap- 
proximation is not violated, namely u^ <C cjpi. Therefore, 
we expect the various tones to reflect that dynamics. 

Figure [TOTa) shows the folded stability diagram, drawn 
in the plane of <l>^^ and <I>^^ using Eq. ([12]) and /3l =45. 
The solid and dashed lines represent stability thresholds 
for positive and negative polarities of <l>^^, respectively. 
Four pairs of solid-dashed lines are drawn, each corre- 
sponds to a different number of flux quanta trapped in 
the SQUID. The black bold line marks the threshold be- 
tween the PNDSZ and the PDSZ. Namely, for excita- 
tion amplitudes |^x^| smaller than the black line, the 
SQUID will reach a LSZ after a flnite number of exci- 
tation periods. On the other hand for |$x^| higher than 
the black line the SQUID will periodically jump between 
LSZs. Note that in this method of excitation, where the 
injected current is kept much smaller than the critical 
current, the SQUID would not be driven to the oscilla- 
tory zone, even by an arbitrarily large ^^^. 

Figure [9fb) shows the simulated second voltage har- 
monics across the SQUID, i.e. at frequency u = 2a;px, 
and for Auj = 0, obtained using E42 parameters. The 
black contours, which are similar to the ones in Fig. 
[TOT a), mark the stability diagram. The distinction be- 
tween the periodic non-dissipative and periodic dissipa- 
tive static zones is clear, and is marked by the white 
bold PDSZ threshold. Like before with the SQUID of 
E38, the PDSZ is characterized by relatively high SQUID 
voltage, which is divided into diamond shaped regions. 
The PNDSZ, unlike with E38, is also characterized by 
strong response along the stability contours. When one 
such contour is crossed, a single transition between LSZ 
occurs. After this transition no additional transitions 
would periodically occur. Nevertheless, these boundaries 
are detectable due to the fact that the inductance of the 
SQUID before and after a transition is different, and also 
due to the fact that the excitation frequency is high, thus 
the effect of changes in the SQUID inductance is measur- 
able. Only transitions across the solid stability contours 
are observed in the colormap of Fig. [9l The reason is 
the sequence of the measurement, which includes sweep- 
ing the DC flux monotonically up and down in the in- 
ner simulation loop. The colormap is obtained while the 
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FIG. 9: (Color online) Parametric excitation experimental 
(a) and simulation (b) results, (a) Reflected power of the 
second voltage harmonics as a function of increasing (blue) 
and decreasing (red) DC flux, (b) Simulated second voltage 
harmonics in the plane of $5^ and $x^. The black contours 
mark the corresponding stability diagram drawn for /3l =45. 
White bold contours mark the threshold to the PDSZ. 



flux amplitude is increased, thus transitions over the solid 
contours are recorded; whereas the decreasing section is 
only used to maintain the consecutiveness of initial con- 
dition, similar to the experiments. 

Figure [9fa) shows measurement results obtained from 
E42. A probe tone having frequency of uj^ is injected into 
the resonator, and the reflected power of the second order 
sideband, i.e. the tone at cJx + Act;, is measured, while the 
DC flux is swept up (blue curve, marked by up headed 
arrow) and down (red curve, marked by down headed 
arrow), and while keeping the AC flux at a flxed ampli- 
tude. Following the blue curve from bottom to top, the 
reflected power experiences a resonance-like absorption 
followed by a saw teeth pattern. The resonance absorp- 
tion pattern originates from the tuning of the resonance 
frequency of the resonator relatively to the frequency of 
the probing signal, thus effectively sweeping the probe 
in and out of resonance. This sweeping is caused by the 
SQUID, which has flux dependent inductance [34], even 
when stuck in a single LSZ. After the SQUID is driven 
across a stability threshold it falls to a new LSZ, into a 
location which is one flux quantum away from the corre- 
sponding stability contour. Thus, if the external flux is 
further increased, it would further drives the SQUID to 
the same direction, and additional transitions would oc- 
cur, spaced apart by one flux quantum. If, on the other 
hand, the direction of the sweep changes (red curve), the 
SQUID would flrst have to be driven across a whole LSZ, 
until reaching the opposite stability threshold. There- 
fore, no matter where the flux sweep changes direction 
along the saw teeth pattern, the reflected power would 
experiences a resonance-like absorption followed by saw 
teeth pattern. 
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Looking back at Fig. [9l^b), the sweeping of ^^^ up 
and down usually drives back the SQUID to its origi- 
nal LSZ. When the AC flux amplitude <l>^^ is increased, 
the borders of a LSZ enclose one another. As a result, 
the sweep across the first stability zone becomes shorter, 
and hence the resonance absorption pattern becomes nar- 
rower along the DC flux axis. This narrowing adds saw 
teeth along the sweep range, and in addition, the SQUID 
might no longer return to its original LSZ, but rather to 
a new LSZ, which correspond to a change of one in the 
number of trapped flux quanta. This creates the sharp 
transition observed in the resonance-like absorption pat- 
terns along the horizontal axis. 



B. Temperature dependant critical current 

When a SQUID is embedded in a resonator, the current 
/x flowing through the SQUID is driven by the resonator, 
and its frequency, thus equals to one of the resonance 
frequencies of the resonator. These first few frequencies 
are only about two to three orders of magnitude lower 
than the plasma frequency of the SQUID [34]. Further- 
more, the heat transfer rate corresponding to hot-spots 
in the NBJJs may be of the order as those frequencies 
or even slower [40] , and therefore effect of heating on the 
dynamics of the SQUID becomes significant. Assuming 
for simplicity that the temperature T/^ (A: = 1, 2) in each 
junction is uniform. The dependence of the critical cur- 
rent on the temperature is given by [53| 



7 — =y[^k) = ^TTVT, 



(13) 



where Icok is the critical current of k^ NBJJ at base 
temperature Tq of the coolant, Qk = Tk/T^ is normal- 
ized temperature of the NBJJ (with respect to its critical 
temperature), 9o = Tq/Tc , and where the function y is 
given by 



2\3/2 



5^(9) = (1-9^)^/^(1 + 9 



^2n1/2 



(14) 



The k (/c = 1, 2) NBJJ heat balance equations read 

Ck^=Qk-Hk{Tk-To), (15) 

where Ck is thermal heat capacity, Qk = V^/Rj is heat- 
ing power, and H^ is heat transfer coefficient. By using 
the notation /Sck = 27rC/eTc/$o^o and Puk = HkCkUJ^, 
Eq. [T5l becomes 



0fe = l^7fe - PHk {Ok - So) ■ 

PCfc 



(16) 



Figure dni panels (b) and (c), show simulation results 
of the SQUID voltage and /3l in the time domain (b), 
and 7+ and 7_ in the phase space (c), calculated for 
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FIG. 10: (Color online) (a) Stability zones in the plan of 
the control parameters $5^ and $x^, drawn for /3l =45. The 
curves are drawn using solid lines for $5^ > and dashed lines 
for $5^ < 0. The numbers beside each pair crossing represent 
the total magnetic flux quanta trapped in the SQUID in the 
corresponding LSZ (in units of $o)- The bold black curve 
marks the PDSZ threshold. Simulation results of the SQUID 
voltage and /3l in the normalized time domain (b), and 7+ 
and 7- in the phase space (c). The simulation is calculated 
using the set of control parameters marked by the plus sign in 
panel (a). Time is normalized by the period of probing tone 
injected to the resonator, Tx, 



the control parameters marked by the plus sign in Fig. 
[TOlfa). The temperature in this simulation is not held 
at the base temperature but rather evolves according to 
Eq. ([16]). The parameters pck = 320 and /^h/c = 5 x lO""^ 
(k = 1,2) were calculated analytically according to Refs. 
[M-ISq. Further explanations about parameters calcula- 
tion is found in appendix B of Ref. [3^ . One expects that 
for the chosen control parameters the SQUID would oscil- 
lates between two LSZs. The dynamics plotted in panel 
(c), however, shows that this is true only during the first 
excitation cycle (dashed line), and that afterwards the 
SQUID oscillates between four LSZs. Such a change in 
the dynamic behavior of the SQUID is also observed in 
the time domain voltage trace plotted in panel (b). The 
number of voltage spikes increases from one during the 
first excitation cycle (dashed line) to three in the end of 
the second excitation cycle (solid line). This behavior 
can be explained by the dynamical change in the value 
of /3l, plotted as red curve in panel (b). The hysteretic 
parameter /3l experiences relaxation oscillations having 
their mean value changing during the first three excita- 
tion cycles. The relaxation oscillations are driven by the 
voltage spikes which dissipate energy and produce heat. 
This heat increases the temperature of the NBJJs, which 
in turn decreases their critical current according to Eq. 
([13]), thus decreasing the value of /3l. The reduction in /3l 
results in a shift of the stability diagram towards the ori- 
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FIG. 11: (Color online) Parametric excitation experimental 
results. The colormaps show the reflected power off the res- 
onator in the plane of $5^ and $x^. Panels (b) and (c) show 
the fitting of the stability diagram to the measured data using 
/3l = 45 and /3l = 35, respectively. 



gin, and consequently the given set of control parameters 
effectively drives the SQUID between increased numbers 
of LSZs. 

Figure [TT] shows the reflected power of the second or- 
der mixing product, i.e. the tone at uj^ + Acj, measured 
in the plane of ^^^ and ^^^ from E42. This measure- 
ment corresponds to the simulation shown in Fig. [9l 
The figure contains three panels, where panels (b) and 
(c) show partial sections of the main colormap shown in 
panel (a), each with a corresponding fit of the stability 
diagram. Looking at panel (a), most of the observations 
agree well with the simulated results, seen in Fig. [9l The 
PDSZ is characterized by strong reflection response and 
stability regions having diamond shapes. The PNDSZ 
experiences strong response along the boundaries of the 
LSZs. The response lines slightly bend for low excita- 
tion amplitudes {^^^ < 3^o), but follow rather strait 
hue for higher amplitudes. Good fitting of the stability 
diagram to these lines is achieved with /3l = 45, as shown 
in panel (b) . The bending of the boundary lines at low 
amplitudes might be due to slow rise of the average tem- 
peratures of the SQUID. The reason why only transitions 
across sohd stability contours are measured is the mea- 
surement protocol, which includes a sweep of the DC flux 
up and down in the inner measurement loop, while mea- 
surements are recorded only during the incremental dura- 
tion of the sweep. In addition to the emerging of the LSZ 
boundaries in the measurement, also the influence of the 
variation of the SQUID inductance within a LSZ on the 
reflected power is observed. Furthermore, the sharp tran- 
sitions along the horizontal axis, corresponding to DC 
flux sweeps that do not drive the SQUID back to its orig- 
inal LSZ, are observed at <^f /<l>o = [1.4,2.75,4.05,4.74], 
in agreement with the simulation results. 

The stability diagram in Fig. [TTTb), drawn for /3l =45, 



predicts that the PDSZ threshold should zigzag along 
the stability contours in the range of ^^^/^q G [7.5,8]. 
The measured colormap, however, shows that the PDSZ 
threshold passes along a strait vertical line starting at 
point (<l>^^ = -6.25^0, ^r = 5.95<^o). A fitting pro- 
cess of the stability diagram to the PDSZ section itself 
produces /3l=35 (see Fig. [TlTc)). This duality in the 
value of /3l can be understood if changes of the SQUID 
temperature are taken into account. The PDSZ thresh- 
old point exactly coincides with one of the LSZ threshold 
contours. The heat generated in single transition across 
that threshold momentarily decreases /3l. An additional 
transition may be triggered provided that the relaxation 
of the first one lasts long enough, which in turn may 
cause further heating of the DC-SQUID. Eventually a 
new mean temperature is achieved for which /3l = 35. 
The measurement protocol dictates that this new tem- 
perature would be kept and that the DC-SQUID would 
stay in the PDSZ for the rest of the measurement. Note 
that the initial value of /3l = 45 differs from the value of 
/3l = 80 that was measured using lock-in amplifier. This 
can be explained by the local heating that the high fre- 
quency flux excitation induces in the DC-SQUID, espe- 
cially in the NBJJs, through the dissipation of circulating 
current due to RF surface resistance. Such reduction in 
/3l from 80 to 45 may be induced by a change of the local 
temperature by approximately 2 K. 



C. The Case of RF SQUID Parametric Excitation 

The stability diagram for a DC-SQUID can be evalu- 
ated numerically, or be analytically approximated for the 
extreme cases of /3l < 1 and /3l > 1. For RF-SQUID, 
on the other hand, it could be exactly evaluated analyt- 
ically. Consider a RF-SQUID having self-inductance L, 
critical current /c, and externally applied magnetic flux 
^x = (^o/27r)^x- The dynamics of the total magnetic 
flux ^ = (<l>o/27r) (j) threading the RF-SQUID loop is gov- 
erned by the potential energy U = UqUjif [57], where 



^RF = (^ - 0x) - 2/3l cos 0, 



(17) 



A local minimum point of i^rf 



and Uo = ^1/ {Stt^L) 

is found by solving S-urf/ dc/) = and requiring that 

d'^u/dcjP' > 0. Clearly, if (p^n is a local minimum point 

of the potential 2Xrf with a given (/)x, then 0^ + 2n7r is 

also a local minimum point of the potential i^rf with an 

externally applied flux of 0x + 2n7r, provided that n is 

integer. 

Lose of stability occurs when d'^u/dc/)'^ = 0, namely 
when cos^ = — 1//3l. This can occur only when /3l > 1, 
since otherwise the system is expected to be monostable. 
This condition is satisfied when (/)x = </>x b, where 



Z>x,b 



TT — arccos 



/?L 



^-y^ 



I - 1- 



(18) 
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FIG. 12: (Color online) Parametric excitation experimental 
results using E19, which has an integrated RF- SQUID. The 
colormap shows the reflected power off the resonator in the 
plane of $5^ and $x*^- The black contours mark the corre- 
sponding stability diagram. They are applied only on half 
of the colormap in order to leave some data uncovered. The 
inset shows a SEM image of a RF SQUID. 



For a general integer n, the local minimum of the po- 
tential itRF near (j) = 2n7r remains stable in the range 



2n7r < (b^ < 



>x,b 



• 2n7r. 



(19) 



Similarly to the case of DC-SQUID, where ^^ is given 
by Eq. (p!Q|) . the boundary contours in the plan of ^^ 
and ^^ are given by the two equations 



^dc ^ ^ac ^ ^ I ^ 



0x,b 

27r 



but are reproducible in measurements. They are only 
detected in measurements with resonators having either 
DC or RF SQUIDs, but not in measurements done di- 
rectly with DC-SQUIDs. Our model only handles the 
DC-SQUID equations of motion, and thus cannot provide 
full description of the dynamics of the resonator-SQUID 
system. Further theoretical work is needed for modeling 
combined TLR- SQUID system, in order to fully under- 
stand these experimental results. 



V. CONCLUSIONS 

In conclusion, we have studied the response of a nano- 
bridge based SQUID embedded in a superconducting 
microwave resonator, nano-bridge based SQUIDs are 
usually characterized by high critical current, and thus 
enhanced metastable and hysteretic response. Several 
phenomena were observed, including periodic dissipative 
static zone in which periodic transitions between local 
stable state occur; hybrid oscillatory zones, in which the 
SQUID is driven to the oscillatory zone by one polarity 
of the excitation amplitude but not for the other, and dy- 
namical variations in /3l due to the effect of self-heating. 
The behaviors of the SQUIDs were compared with theory 
both analytically and numerically with good agreement. 
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for the largest and smallest values of ^x, respectively. 

Figure [12] shows the reflected power of the third order 
mixing product, i.e. the tone at uo^ + 2Aco', measured in 
the plane of <l>^^ and ^^ with El 9, which has an embed- 
ded RF-SQUID instead of a DC-SQUID. The black con- 
tours represent the stability diagram, plotted using Eqs. 
(|2Q|) for /3l = 1.5. The stability diagram qualitatively 
matches the measured data for the first few stability di- 
amonds. 

Note that the value of /3l for the RF-SQUID of E19 is 
more than order of magnitude smaller than that of E38 
and E42. This might be due to the fact that this RF- 
SQUID is fully fabricated on Silicon Nitride membrane, 
and thus has reduced thermal coupling to the coolant 
compared to the DC-SQUIDs of E38 and E42. This in 
turn might lead to an increased local temperature of the 
NBJJ, and to a degraded value of /3l. 

Note also that the diamonds themselves have distinct 
reflection patterns inside their bounded zone (see Figs. 
[TTI and [T2]) . which are not perfectly periodic with ^: 



Appendix A: Nano-bridge current-phase relation 

The CPR of a single short channel of transmission r is 
given by ^] 



-t^<^'' 



where 



j{i) 



r sm7 



1 — r sin^ ^ 



(Al) 



(A2) 



ac 

X •) 



The NBJJs in our devices are not ideal one dimensional 
point contacts as Ref. [50] assumes, however, we have 
found that the above simple analytical result resembles 
the CPR which is obtained by solving the Ginzburg- 
Landau equation in the limit of short bridge (in com- 
parison with the coherence length) [24]. Let 70 be the 
point at which the factor J (7) has its largest value J (70), 
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which is given by 



J (70) = 2 Y^2 (1 - VT^) - T . (A3) 

Using this result the current / can be written in terms of 
/c as I/Ic = F (7), where 



and 



72 + /^d72 + (1 - (^0) y (©2) F (72) 

--S— (71 - 72 + 27r^x/^0) = I /ho + ^n2. 

Plo 



(A6) 



F(7) = 



r sin7 



2 J2 (1 - VT^^) - tJi - Tsin^ § 



• (A4) 



Replacing the sin 7k terms {k = 1,2) in Eqs. Q and (|8j) 
by F (7k) leads to the following modified EOMs 



71 + /^D7i + (1 + c^o) ^ (©1) F (71) 

+ ^^ (71 - 72 + 27r<l>x/<^o) = /x//cO + ^nl, 
PLO 



(A5) 
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